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Abstract 

Investigations of two resonant planets orbiting a star or two resonant 
satellites orbiting a planet often rely on a few resonant and secular terms 
in order to obtain a representative quantitative description of the system's 
dynamical evolution. We present a semianalytic model which traces the 
orbital evolution of any two resonant bodies in a first- through fourth-order 
eccentricity or inclination-based resonance dominated by the resonant and 
secular arguments of the user's choosing. By considering the variation of 
libration width with different orbital parameters, we identify regions of 
phase space which give rise to different resonant "depths," and propose 
methods to model libration profiles. We apply the model to the GJ 876 ex- 
trasolar planetary system, quantify the relative importance of the relevant 
resonant and secular contributions, and thereby assess the goodness of the 
common approximation of representing the system by just the presumably 
dominant terms. We highlight the danger in using "order" as the metric 
for accuracy in the orbital solution by revealing the unnatural libration 
centers produced by the second-order, but not first-order, solution, and by 
demonstrating that the true orbital solution lies somewhere "in-between" 
the third- and fourth-order solutions. We also present formulas used to in- 
corporate perturbations from central-body oblateness and precession, and 
a protoplanetary or protosatellite thin disk with gaps, into a resonant sys- 
tem. We quantify the contributions of these perturbations into the GJ 876 
system, and thereby highlight the conditions which must exist for multi- 
planet exosystems to be significantly influenced by such factors. We find 
that massive enough disks may convert resonant libration into circulation; 
such disk-induced signatures may provide constraints for future studies of 
exoplanet systems. 

Keywords: RESONANCES, ORBITS, PLANETARY DYNAMICS, EXTRASOLAR PLANETS, PROTOPLANE- 
TARY DISKS 
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1 Introduction 



Resonances have played an increasingly important role in the dynamical analysis of planetary 
and satellite systems. The ongoing discoveries of extrasolar planets, Solar System satellites, Kuiper 
Belt Objects, and asteroids have sparked a resurgence of interest in resonant systems. The variety 
of resonances observed or thought to exist in nature showcases the utility of a versatile model which 
can help determine what approximations are sufficient or inadequate for future detailed studies. 

The term "resonance" typically refers to any system which features a commensurability of fre- 



quencies. 
resonance ( 



hese frequencies may refer to an object 's electromagnetic forcing, the so-called "Lorentz" 



Burns et al. 



Burns et al 



1985 



Schaffer and Burns. 



2004), a n object's spin ( 



Biasco and Chierchia 



1965a 



Peak 



1976 



2002 



1986 



1992 



Hamilton and Burns 



Goldreich and Peak 



Flvnn and Saha 



Greenberg 



1977 



1966 



Murdock 



2005), or an object's orbita 



Malhotra 



19941 : 



Morbidelli 



1993: 


Hamilton 


. 1994 


1978 


Celletti. 


1993 


elements f 


Goldreich . 



20021 ). Commensurabilities 



between orbital frequencies of any number of bodies may exist. Most known cases of orbital reso- 
nance occur between two objects revolving around a massive central object; our model evolves such 
two-body "orbit-orbit" resonances. Two-body orbit-orbit resonances have been subdivided into a va- 
riety of designations, but naturally are split into "mean-motion" and "secular" resonances. Secular 
resonances occasionally occur inside mean-motion resonances, but not vice-versa. "Secondary" reso- 
nances, which arise when the libration frequency of the primary resonance is commensurate with the 
circulation frequencies of the secular mode of motion, also lie just inside mean-motion resonances. 
Some secular resonances include "eccentricity-type" resonances, which involve commensurabilities 
among longitudes of pericenters, "inclination-type" resonances, which involve commensurabilities 
among longitud es of ascending nodes, and "m ixed" resonances, which include relations among both 

1999 . p. 359)[. 



types of angles (Murray and Dermott 



No orbit-orbit resonances exist among the 8 planets of the Solar System, in contrast to the 
several resonances present in recently discovered multiple-planet extrasolar systems. Jupiter and 
Saturn, however, are in a 5:2 near-resonance known as "The Great Inequality", which has been 



known since the time of L aplace and recently studied in 



Michtchenko and Ferraz-Mellol (120011 ) and 



Franklin and Soperl (120031 ) . Relative to known exosystem resonances, Solar System resonances 



Wakhidovl l|200lh classifies three-body resonances among Jupiter, Saturn and an asteroid as "mixed." 
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involve comparatively small eccentricities, are better identified due to their proximity to Earth, and 
have undergone scrutiny for a longer period of time. The ongoing discovery of extrasolar planets 
has led to at least 20 confirme d multiple-planet exosystemsy almost half of which (see Figs. 1 and 



2 of 



Veras and Armitage 



20071 ) reside close to a lst-4th order mean motion resonance. GJ 876 b 
and c are two of the most well-studied resonant exoplanets, and are widely thought to reside in the 
str ong 2:1 resonance. 



Marcy et al. 



(1200 ll ) reported the presence of two planets in the GJ 876 system with orbital 



parameters suggestive of the presence of a mean motion resonance. The system has since become a 



catalyst for more careful studies of the 2:1 resonance, w^ 



for stability and possible prior evolution (IBeauge et al. 



r ich re pr esents a crucia l 



2006). 



geometry of the 2:1 resonance as app lied to GJ 876, an c 



Psychoyos and Hadjidemetrioul (120051 ) and 



Beauge et al. 



Ferraz-Mello et al. 



dynam ical marker 



Lee and Peald (120021) exp l ore th e 



(I2003h 



Lee! (J200J) 



(120061 ) sample the resonant phase space 



further in order to demonstrate the diversity of the resonant configurations. Orbital fits suggest 
that the dominant resonant and secular angles all librate about 0° with well-defi ned amplitudes, 



and a ll ow for both res onant planets to harbor a modest nonzero mutual inclination (jLaughlin et al. 



2005|). 



Ji et al. 



(120021 ) perform 1 Myr simulations of the GJ 876 system, assuming different initial 



relative inclinations. 



Ji et al 



(120031 ) suggest that the planets in GJ 876 are likely to undergo apsidal alignment and 
derive a criterion f or determining if a parti cular secular argument in a two-planet system is in libra- 



tion or circulation. 



Snellgrove et al. 



(120011) and 



Kley et al 



( 120051 ) use hydrodynamic simulations to 



explore the possibility that the currently observed eccentricities of the GJ 87 6 pla nets are the result 
of differential migration from the nascent protoplanetary disk. 



Ward 


(1981 


) and 


Kiev et al. 


(2005) 



give fo rmulas for the apsidally induced precession of a planet due to a disk. Although 



Rivera et al. 



(120051 ) suggest the presence of a third planet in the GJ 876 system orbiting at ~ 0.02 AU, the 
planet's small (0.23mj up ) mass and circular orbit are unlikely to affect significantly the resonance 
between the other planets. 

A salient feature of orbit-orbit resonances is the domination of just one or a few dynamical terms, 
often dubbed "resonant terms" or as "secular terms" (depending on the nature of the terms), in 
the system evolution. Resonant models for particular systems are often built around the terms 



2 From the on-line Extrasolar Planets Encyclopedia, at http://vo.obspm.fr/exoplanetes/encyclo/catalog.php 



chosen. Our model can help investigators determine which terms dominate a particular system by 
allowing the user to include up to twenty resonant and/or secular ar guments for each simu l ation 



Sometimes, however, additional effects, such central 



p. 264-270) or central-body precession (iRubincam 



- body oblateness ([Murray and Dermott 



1999 . 



2000l ). can play a crucial role in the system 



evolution. Further, in the formation stages of a planetary or satellite system, a disk may be 
present. The mass and breadth of the disk then can strongly influence the resonant dynamics. Our 
code thus includes the ability to toggle the effects of oblateness, precession and a nascent disk. 
Our approach to evolving the resonant bodies involves retaining a classical set of orbital elements 
while simultaneously tracking each resonant angle and its time derivative. The method allows for 
additional effects to be incorporated through the system's potential. By treating the dynamical 
equations in as general manner as possible, our code is effective for a wide range of initial orbital 
parameters and masses, but cannot accurately model crossing orbits, high eccentricities, nor high 
inclinations. 

In Section El we derive our core model, connect its analytical and numerical aspects, and relate 
our formulations to action-angle variables and system constants. Section [3] begins the resonant anal- 
ysis by describing asteroidal motion in the restricted three-body problem, and Section H] illustrates 
how libration width varies with a variety of parameters for two massive planets. We then apply 
the model to a real exoplanetary system, GJ 876, in Section [5J with a term-based perturbative 
treatment. In Sections [6l [7] and [U we present both averaged and unaveraged formulas expressing 
the gravitational effects of central-body oblateness, precession, and a thin disk, and detail their 
contributions to GJ 876, thereby illustrating the necessary conditions for their consideration in 
other exosystems. We conclude in Section [91 

2 The Core Resonant Model 
2.1 Definitions 

The orbital motions of both resonant bodies in an isolated system a re described completely by 



"Lagrange's planetary equations", which (jBrouwer and Clemence 



1961 



p. 273-284) derive without 



approximation. Despite their name, these equations are not restricted to describing planetary 
motion, and are dependent on two "disturbing functions," which represent potenti als that arise from 



applying Newton's gravitational force laws to the central and resonant bodies. 



Ellis and Murray 



(120001 ) summarize the history of the development of the disturbing function, and help describe 
its modern usage. The disturbing functio ns are infinite linear combinations of cosine terms with 
arguments of the form (IKaulal . Il96ll . Il962l ). 



(f)(t) = j 1 \ 1 (t) + j 2 X 2 (t) +jl, zu TU 1 (t) + j 2> ^ZJ 2 (t) +ji,n^i(t) + j 2 ,Q^ 2 (t), 



(2.1) 



where A represents mean longitude, w represents longitude of pericenter, Q represents longitude 
of ascending node, the "j" values represent integer constants, and the subscript "1" refers to the 
outer planet while the subscript "2" refers to the inner planet. Secular arguments are those for 
which ji = j 2 = 0, and arguments which are multiples of one another are distinguished because 
they have different coefficients whose magnitudes may vary drastically (by orders of magnitude). 
The set of elements which define each resonant argument are subject to constraints known as the 
d'Al embert rela t ions, which require th at j\ + j 2 + ji ro + j 2w + jxu + j 2 ,n = and j\ n + 32 a is even 



(e.g. 



Greenberg 



1977 



Hamilton 



19941 ) . 



Typically, orbital elements in celestial mechanics are, by default, assumed to be osculating, i.e., to 
satisfy the Lagrange gauge (or Lagrange constraint), which demands that the elements parameterize 
instantaneous conies tangent to the perturbed orbit. Under this constraint, the dependence of 
the perturbed velocity upon the elements has the same functional form as the dependence of the 
unperturbed velocity upon the elements. Derivation of the standard planetary equations in the 
forms of Lagrange or Delaunay is based on this constraint. This derivation, however, also exploits 
the assumption that the perturbations depend solely upon positions, not upon velocities. In the case 
of velocity-dependent disturbances, the planetary equations for osculating elements assume a more 
complicated form. Specifically, they acquire new terms that are not parts of the disturbing function. 
For the reason, osculating elements become very mathematically inconvenient when considering 
velocity-dependent perturbations. Perturbations of this type arise in problems with relativistic 
corrections or with atmospheric drag. They also emerge when we switch from an inertial reference 
frame to a precessing one - below we shall encounter exactly this situation. 

Fortunately, even under velocity-dependent perturbations, we can restore the standard form of 



the planetary eq uations. This, however, c a n be achiev ed on 



demonstrated by 



Efroimsky and Goldreichl ( 2003 



20041 ) and 



y by sacri fi cing th e osculation. As 



Efroimskyj (j2005al Jbl). there exists a 
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non-Lagrange gauge which returns to the planetary equations their customary form. The advan- 
tage of this approach is that even the velocity-dependent perturbations appear in these equations 
simply as parts of the disturbing function. The consequence is that the elements rendered by these 
equations are no longer osculating. Such elements model the orbit with a sequence of conies that are 
not tangent to this orbit. Thereby, these elements return the right position of the satellite but not 
its correct velocity. They are called "contact elements" (term offered by Victor Brumberg) . Though 
these elements were rigorously defined and comprehensively studied only very recently, they had ap- 
peared in the thitherto literature whenever someone incorporated a velocity-dependent perturbation 
into the disturbing function and then substituted this function into the standard planetary equa- 
tions. By performing this sequence of operations, one tacitly postulated a certain non-Lagrange 
gauge, i.e ., accepted a cert ain "amount of nonosculation . " Fo r the first time, this s ituation oc- 



curred in 



Goldreichl (Il965bl ) and later in 



Brumberg et al 



(119701) and 



Kinoshital fll993[ ). Goldreich 



and Brumberg noticed that the elements furnished by these written equations were nonosculating. 

Though the contact elements differ from the osculating ones already in the first order (over 
the perturbation caused by the transition to a precessing reference frame), the secular parts of 
contact elements differ from those o f their osculating counterparts only in the second order. This 



result, proven by lEfroimskyl (l2005bl ). is valid only in the case of uniform precession and only for a 
solitary satellite, in the absence of any other disturbances (like the gravitational pull of the Sun or 
of another satellite). However, numerical simulation has shown that even in realistic situations of 
variable precession, the deviations between the secular parts of the contact and the cor responding 



osculating elements accumulate very slowly, for a solitary satellite (jGurfil et al. 



20061 ). Thus, in 



practical cases, one may safely assume that for a solitary satellite, the secular parts of the contact 
elements make a very good approximation to those of the appropriate osculating ones. For our 
model, the velocity correction between both sets of elements is on the order of the relativistic 
correction, an effect orders of magnitude smaller than any we consider. 

We define the "order" of a resonant argument as \ji + j 2 \, and denote each resonant argu- 
ment by the set {ji, j 2 , j 2 ,zu, ji,n, J2,n}- The mean longitude is directly proportional to mean 
longitude at epoch, denoted by e^, such that \k = zuk + Mk = vr^ + fi^ + Mk + f* Q nk(t')dt' = 
e k + J* [i^ 2 ^ a k(t'Y~ 3 ^ 2 " > dt' , where M k denotes mean anomaly, 7r fe denotes argument of pericenter, 
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and n k denotes mean motion, with k — 1,2. Henceforth, the letter "/c" will represent the dummy 
variable which denotes the outer and inner resonant bodies. The mean motion is related to its 
semimajor axis and mass through Kepler's third law by n\a\ = fi k , where fi k = Q (mo + rn k ), with 
m k representing the planet's mass, mo the central body's mass, and Q the universal gravitational 
constant. Using the above result in conjunction with the time derivative of Equation (12.11) (all time 
derivatives will henceforth be denoted as overdots) yields: 



4>(t) = h$<h 5 (t)+jih(t)+j2^a 2 2 (t)+j2e 2 (t)+j 1 ^w 1 (t)+j2^w 2 (t)+j 1) nQ 1 (t)+j2,nn 2 {t). (2.2) 



Equation (12. 2p provides a crucial component of our model; we will later show (Eq. 12.161) that the 
right-hand-side may be expressed in terms of <ft, a^, the eccentricity e k and the inclination I k only. 



2.2 Equations of Motions 



The complete set of Lagrange 's planetary equations, with t he disturbing functions denoted by 



1Z k , can be manipulated to read (IBrouwer and Clemence 



1961 



p. 284-286), 



where, 



da k i &R k 

-dT = akAk > 1 ^ (2 ' 3a) 

It = Ak * Ak ^ k ~ ^ (2 ' 3b) 

dl k -i 8TZ k -i &R k -i on k 

~TT — ~ a k A kt 4;A k:5 — a k A k> ±A k £— a k A k ^A k $ , (2.3c) 

dt o\ k Ow k o\i k 

de k - &R k —- dK k ~- dTZ k 

— = -al A k l — h a k 2 A k)2 A ki3 — h a k 2 A kA A k $— — , (2.3d) 

dt oa k oe k ol k 



dw k _ -i dTZ k -i dTZ k 

— — — a k A k 2 ~z V a k A k AA k 5— : 

dt de k K dl, 



a k 2 A kt2 ^— + a k 2 A kA A k)5 -^- , (2.3e) 



dQu - 1 - dTZi, 

■W^'A^Hit' (2 ' 3f) 

Afc = / n k dt + e k , (2.3g) 
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A k ,i 


= 2^ 


2 

? 


A k>2 


= A** i 




A k ,3 


= 1 - 




A k ,4, 


= A*»* 






= tan 






= CSC 


4). 



(2.4) 



A primary feature of our model is the ability to include whichever, and as many, terms in the 
disturbing function that one wishes depending; on the resonant si tuation. Our heliocentric disturbing 



( H) 

function, lZ k , may be written as (Murray and Dermott 



1999 



p. 329): 



r'E 

i=l 



Ui 



p=l 



cos< 



(0 



(2.5) 



such that i labels each argument 0, and Ui — 1,3 or 11 because the expansion is taken to fourth- 
order. We represent the {0,0,0,0,0,0} secular argument, for which U = 11, as K, and any 
ot her secular or res o nant argument as h. The expansion of the disturbing function is that due 



to 



Ellis and Murray 



J2000|). 



The quantity X^'^ represents the infinite linear combination of powers of eccentricities and semi- 



inclinations that accompany each term, and the quantities C^' p ^ are functions of the semimajor axes 
alone. The summation in the brackets indicates the possible presence of more than one coefficient 
associated with a single resonant or secular argument. 

The form of the disturbing function in Eq. (12. 5p fails to model accurately resonant bodies 
in crossing orbits due to the singularity at orbit inter section, and the resulting se ries convergence 



1999 



p. 250). There- 



depends on the orbital proximity to this crossing point (IMurray and Dermottl . 
fore, this model cannot reproduce the evolution of resonantly locked orbit-crossing pairs such as 
Neptune and Pluto. Further, the convergence domain of the expansion on which Eq. (12.51) re- 
lies precludes realistic soluti ons for high eccentricities. Although high e ccent ricity expansions of 



the distu rbing function exis t 



utility of 



(IRoig et al. 



1998 



Beauge and Michtchenko 



20031 ). we investigate the 



Ellis and Murray! (120001 ) traditional expansion about zero eccentricities and inclinations. 



The expansion converges only for e < 0.66; however, the Sundman criterion, applied in Section [5], 
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restricts the magnitude of the eccentricities even more. We assume, 



where for p > 1 , 



/. ,x li'W I h' W I h'W I h' w 

vY — 61 611 



(2.6) 



92 
2 ■ 



"p— 1 



01 
02 



(p + 3) mod 5 
2 

(11 — p) mod 5 



where the "semi-inclinations" are Si = sin(/ 1 /2) and S2 = sin(/ 2 /2). 
functions of the masses and semimajor axes alone: 



[2.7) 



X8) 



The quantities are 



u 2 



Qm 2 



-2 , 

" Jint > id 



(2.9a) 



^2.9b) 



where a = 02/01, and /g*f ; are constant "indirect" "internal" and "external" contributions, 
and f^' p ^ is the "direct" contribution, which is a function of ji and a. The indirect contributions are 
zero for the majority of orbital resonances. In many resonant stu dies, f^'^ is treated as constant by 
setting a as a constant computed from the initial semimajor axes; Ferraz-Mello ( 1986)1 ) demonstrated 



the danger in doing so, and hence we do not make that assumption here. We use 



t p) = ff p) Ht)) = Y,^ p) *(t) 1 

1=1 



(2.10) 



where ftf' p ^ are constants specific to each term (i,p), and are what give each resonance its unique 
character. The Appendix describes our method for obtaining k[ 1 ' p ^ values. Our code incorporates 



K 



ihP) 



coefficients for the and C^' 11 ' terms corresponding to all unmixed first thru fourth-order 



resonances, and for (zeroth-order) secular terms up to degree 4 in eccentricities and inclinations. The 



10 



disturbing function in Eg. (12.51) assumes helio c entric coordinates, such that 1Z X /TZ 2 H ^ = m 2 /mi. 



In Jacobi coordinates rtBrouwer and Clemence 



1961 



p. 589), 



7l[ J ^ mQm2 (mo + mi + m 2 ) 



K 



(■') 



m x (mo + m 2 y 



(2.11) 



By using Eq. (12.51) and Ck = cos (Ik/2), we can take partial derivatives of the disturbing functions 
and insert them into Lagrange's planetary equations, which yields: 



dak 
~~dT 

de k 
dt 

dlk 
dt 

dei 
~dt 



de 2 
~dt 

dwk 
~dT 

dQk 
~dT 



r'E 

i=l 
l 01 

i=l 
i 00 

i=l 
l 00 

i=l 

oo 

— - X — ■> 

i 00 

i=l 

1 °° 

i=i 

i °° 



L k B k,l 

p=l 

Ui 

P =i 
p=l 



sin 1 



sin 



sin 1 



dai 



i=l 



E 

P =i 

Ui 

P =i 

u i ap{i,p) R (i,p) 

E Q4 -"2,1 

p=l 1 J2 

E/7(*jP)_g(*.p) 

,p=i 
,p=l 



(*>p) 



COS ( 



Jl 



(0 



COS ( 



cos<x> v ; + a 9 a 



r'E 



i=l Lp=i 



(i,p) r(«>p) 

2,4 



COS 



COS 



cos< 



(2.12a) 
(2.12b) 
(2.12c) 
(2.12d) 
(2.12e) 
(2.12f) 
(2.12g) 
(2.12h) 



where 
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r>(i,P) 
D k,l 




(2.13a) 


r(*>p) 
D k,2 


= ^, 2 ^3j|^ (j ' p) +^,2j|;U (i ' p) 5 


(2.13b) 


r(*>p) 
n k,3 


= A^^x^) +igUMA l5 ^ (i,p) +j|U,44^ M . 


(2.13c) 


r(*>p) 




(2.13d) 


r(*>p) 


= ^^e^*) + \hlfA kA A k ^c k X^\ 


(2.13e) 


r(*.p) 


= \htfA kA A kfi s- k l c k X^\ 


(2.13f) 



with 



,(*.p) 

4 Jfc,e 



and 



•(») 

3 k,vc 



•(*) 
•7 fc,tc 



0k(p)> 



,0, 



+ 2 (l-|fc-p+ 1|) |sgn (btil + ligil) 



(*,p) 



Ll*,P 



•(*) 

.7fc.fi 



.(i) 



0fc(p), 



+ 2(l-|A;-p+l|)|sgn(|j«| + |j 



/2,n 



x and p = 1 
x and p > 1 
N and p > 1 
K and p = 1 , 

x and p = 1 
, x and p > 1 
N and p > 1 



(2.14) 



(2.15) 



k 0, K and p = 1. 

Both the A and -B auxiliary variables are functions of only the eccentricities and inclinations. 

Now combining Eqs. (I2.2p and (2.11d)- fl2.12hj) provides the following compact resonant equation: 

for u = 1, ...Z, where Z denotes the number of cosine arguments retained in the disturbing function, 



(«) = ( D(l,u) cos0 w ) 



(2.16) 



where 



12 



D (i,u) = g 



and, 



P =i 

2 



,-1 / n ~2p(hP) ( ' •(«) r(*.P) I „•(«) d(«,p) I „•(«) r(«,pA „•(«) /l V-(*>p), 
l i Z^\ a fc C fc ^* + 3k,z* B k,5 +3k,ii B kfi )-Jk A k,lX" >i 

k=l 



(*>p) 



<9a fc 



(2.17) 



1 3 

2 



(2.18) 



fc=i 



Equation (12.161) illustrates that for a disturbing function with Z cosine terms, the time derivative 
of each argument is a linear combination of the cosines of all Z arguments, with coefficients that 
are functions of m , mi, m 2 , ai, a,<n e i> e 2> h an d -^2 only. The arguments may be resonant or secular. 
Equations (12.121) and (12.161) represent a self-consistent set of 12 + Z first order coupled differential 
equations, which our model integrates directly. We use an adaptive-timestep fourth-order Runge- 
Kutta integrator with user-defined accuracy parameters. 

2.3 System Constants 

Although some quantities, such as energy and angular momentum, represent constants of any 
isolated physical system, they strictly no longer represent constants when a finite number of terms in 
the disturbing function are used to model the evolution. Further, in the restricted problem, energy 
and angular momentum no longer represent constants of the motion. Regardless, the variation of 
these quantities may be negligible depending on the system and the terms chosen. Further, if only 
one or a few terms are considered, then additional constants may arise, as we detail in this section. 
These additional constants prove useful in studies of, for example, two resonant bodies in which one 
is much less massive than the other, such as the case with a main belt asteroid and Jupiter, Titan 
and Hyperion, or a terrestrial extrasolar planet and a giant extrasolar planet. The system angular 
momentum, c, may be expressed as: 



c i + c i + 2ciC 2 [sin Ii sin I 2 cos (Q\ — fl 2 ) + cos h cos h 



(2.19) 



where, 
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Cfe 

and in Jacobi coordinates, 



= m k Jg (m + rnjfc) a fe (1 - e|), 



(2.20) 



Cl 
^2 



(to + TO 2 ) TOi 

y/m + m 1 + m 2 
m m 2 



Ga x (1 - e\) 



(2.21) 



£ci2 (1 



^m + m 2 

A time-independent Hamiltonian provides the energy constant for the system, and allows one 
to construct canonical sets of variables which can reveal additional constants and provide further 
insights into the considered system. In terms of orbital elements, the Heliocentric Hamiltonian may 
be approximated as: 



n 



H 



Qvtlqvh\ Qm§m 2 m,\ 



oo Ui 



2ai 



2a 2 



a i 



j2J2\ c i' p)xii,p) 



cos- 



(2.22) 



i=i p=i 

One may attempt to avoid the approximate nature of the Heliocentric Hamiltonian by expressing 
the Hamiltonian in Jacobi coordinates (TCj)- The Hamiltoni an for a three-body system in Jacob i 



coordinates, a form of which has been used by several authors 

~1 



1983 



Konacki et al. 



200( 



For d et all 200 



by others (iMalhotra and Dermott 



19901 : 



0; 



Lee and Peald. 



Ferrer and Osacarl . 



Harrington 



1968 



1969 



Sidlichovsky 



2003) and explained from first principles 



19941 ) may be expressed as: 



where 



with 



K 



Qm m 2 G (to + to 2 ) toi 



2ao 



Qm^m\m 2 



2ai 



-K, 



1=2 V r l/ 



m 



i-i 



-m 2 



< 1, 



(2.23) 



(2.24) 



(2.25) 



(to + m 2 ) 

where tp is the angle subtending ri and r 2 . Note that no indirect terms appear in Eq. (12.231) . 
thereby eliminating the need for "internal" and "external" nomenclature. We seek to express Jtj 
in terms of m^, a^, e^, 1^ and (j)^' only. 



As partially demonstrated by 



Lingi (Il99ll ). the use of Jacobi coordinates does no t alter the 



derivation of the disturbing function given by Eq. (6.36) of iMurray and Dermott 



(1999 



p. 232), as 
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the coordinate system transformation does not alter the angle between the radius vectors nor the 
expansions performed on the orbital elements of the indi vidual orbits. Therefore, Eq. (12.241) may 
be expressed as Eq. (6.36) of iMurray and Dermottl (119991 . p. 232), with a replaced by o^, where 



a l M = a l Wli, I > 2. 



(2.26) 



All other variables in that equation remain unaffected. The Jacobi Hamiltonian thus reads: 



Qm m 2 G (rn + m 2 ) m x ^mo™i™2 



2a, 



N Ui 

EE 

i=l p=l 



. 1=2 



cos0 w . (2.27) 



2a\ a\ 

We can establish a set of canonical action angle variables from the Hamiltonian; we choose the 
following set: 



Afc = M k + ir k + Q k 

Ik = -TTfc - 

z k = — Q k 



A k = Mky/Ok' 



(2.2* 



r fe = A fe (l-,/l-ei 



Z k = A k J (1 - e%) (1 - cos 4 



where A4 k is a function of the masses. In the literature, when M. k — VGmnmk/ \/mo 



action-angle variables have been classified as " Poincare variables" (jMurrav and Dermott 



60) "mass-weighed Poinc are elliptic variabl es" (IMichtchenko and Ferraz-Mello 



fied Delaunay variable s" (jMorbidelli 



1976 



Morbidelli 



2002L p 



mfe, the 



19991 . p. 



20011 ) . and "modi- 



■ P- 33), 



.2002 
\ 

•Mk — m kVG m o ( jVaradi et al. 



Mk = VQm (IMorbidelli 



35). In some cases, M. k = \JQ (mn + m k ) 



2001 



Eui Chang and Marsden 



Peald . 



a, 



2003 , or 



19991 ). The form of M. k is chosen based on the Hamiltonian of the 



system. When the Hamiltonian is written in Ja cobi coordinates, M-? = \/Gmnm?J\/ mn + m% 



1968 



Sidlichovsky 



19831 ). With 



while M.\ = \[Q (m + m 2 ) mi/y/m + mj + m 2 (IHarrington . 
7r k + Q k = G7 fe , the angles correspond to those seen in the resonant arguments. A canonical 
transformation to the variables (8i,Qi), I = 1...6, can now be applied to the Hamiltonian such 
that: 
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with 



j?X 1 + j®\ 2 - j« , 7l + [j-> + $> + , , 72 



;(0 _i_ „•« , „•« 



7 i — 



•(*) , •(*) 



#3 = 

o 6 = 



•(*) , •(») ' 
Jl + J2 

Jl M + 32 A 2 - 3 3 -k n ,n Zl + Ul + -?2 + J3-fc n ,Q ) ^2 



Ai 



3i +32 



Jl + J2 



6i = 



e 2 



e. 



3 fc^; 



•(*) 



•(*) 
3 ,vj 



3 fc-j^7 I fc-^j J 



Jl + 32 ) Z 3-fe n 



e 4 = 



.?3-fcfi,fi 
■71 T,7 2 ^J3-fc n ,Q 

•(*) 



e 5 = (j? + 3f) Ai - jf ) (r fcro + r 3 _ fcro + + z 3 _ fen ) , 
©e = (3? + 3f) a 2 - 3? (r fcro + r 3 _ fcro + + z 3 _ kn ) , 



I 7 W 







A: 



'n 



1, 3\ 







(2.41) 



2, otherwise 12, otherwise 

For mixed resonances, A; ro and k^ may be set to 2 when = and j[^ n = 0; the choice is arbitrary. 

I chose this transformation both for application to each Hamiltonian and in order to take into 
account several types of resonances, including unmixed and mixed eccentricity and inclination 
resonances to arbitrary order. The transformation is not well suited for the rare coupled eccentricity 
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and inclination resonances, which, according to the d'Alembert rules (iMorbidelfil . 120021 . p. 35-36), 
must be at least of order 3. For disturbing functions with a single resonant argument, that argument 
should equal a multiple of 8\ or # 3 . For disturbing functions with several resonant arguments, each 
argument may be equated with an appropriate linear combination of the Q\ valueqj. The number 
of possible resonant arguments associated with given values of and is \J® + + lj. 
Typically, when multiple resonant arguments are present in a disturbing function, they all have the 
same values of and ■ This situation allows one to explore the effects of "resonance splitting" , 
when ji and j^ 1 remain fixed as % changes but the other j values change as i does. If j± and 
are fixed for several i, then the canonical transformation chosen above allows for all such arguments 
to be represented as a linear combination of only 81 and 82 or only 83 and 84. 

This transformation allows us to obtain constants of the motion, as all quantities in the Hamil- 
tonian other than the resonant arguments may be expressed in terms of the actions only. These 
constants may be determined immediately from the initial conditions, and then used in the subse- 
quent analysis. When the disturbing function contains a single resonant argument, the Hamiltonian 
does as well, and several constants of the motion exist regardless of the form of Aik- For any 
eccentricity-type of resonance, two of these constants, iV\ and N 2 , may be expressed in terms of 
a k (t) and e k (t) only: 



N k = y/ako~ 



.(i) 

3 fc.tc 



1 



Jk 



(i) 



(2.42) 



where at and ek represent initial, known values. For any inclination type of resonance, the two 
constants may be expressed in terms of Ofc(i) and Skit) only: for k — 1, 2, 




(2.43) 



where Sk similarly represents an initial known value. The constants in Eqs. (12. 42ft and ( 12.431) may 
be derived from either a Hamiltonian or a non-Hamiltonian approach. For the former, all G; except 

3 See, for example, Michtchenko and Ferraz-Mello (2001), who include four resonant arguments in their Hamiltonian. 
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either 9i or O3 must be constants: for eccentricity resonances, both actions are constants, 
while for inclination resonances, both actions are constants. The resulting 5 constants may be 
manipulated to yield Eqs. (I2.42p - fl2.43p . Because T k = Y k {a k ,ek) whereas Z k = Z k (a k , e k , s k ), no 
inclination terms, constant or otherwise, appear in Eq. (I2.43p . whereas constant eccentricity terms 
do appear in Eq. (12.431) . In the non-Hamiltonian derivation, dividing Eq. (I2.12al) by Eq. (12.12bl) . 
separating variables, and integrating will yield Eq. (12.431) . One follows the same procedure with 
Eq. (I2.12al) and Eq. (I2.12cl) in order to yield Eq. (I2.43p . with one difference: the eccentricity in 
Eq. (I2.12cl) must first be written in terms of a k using Eq. (I2.42p . 

The quantities N\ and N 2 are constants in both the heliocentric and Jacobi coordinate systems, 
and are independent of any masses. Each provides a relation among the orbital elements of a single 
body, but does not correlate the parameters of both resonant bodies to one another. A constant 
which achieves such a correlation would differ, albeit slightly, depen ding on the coord inate system 



used, and the form of M. k . If we take M. k = \/Qm m k / \fm + m k ( IMorbidellil . 



2002 



, p. 35), then 



using only the constants 9s and 96 immediately gives the following relation: 



AT _ I J I " U Z 1 

iV 3 = V a 2o _ 771 V a io 



-3i 


m 2 


•( 


;) 


J2 


mi 


-Ji 


m 2 


•0 
J2 


nil 



(2.44) 



a 2 -7TT v a i- 



Whereas for the Jacobi Hamiltonian, 



N 3 = I — — I : 7i 

(m + m 2 ) 2 







m 2 


u 


''mi; 






\32 


mi J 



(2.45) 



02 I -?7i — I ; -3 v a i- 

(m + m 2 y 



Note that in the limit mo ^> nil, m 2 , both expressions for A^3 tend toward equivalence, as expected. 
Also consider that for k = 1,2, in the limit m k — > 0, 03.^ is a constant. Physically, this result is 
sensible, as a negligible mass should not perturb a larger mass off of its original orbit. In the case 
of multiple resonant arguments but no secular arguments, the Hamiltonian can be expressed as an 
appropriate linear combination of either 9i and Q 2 , or 93 and 94. Hence, the other four actions 
are constants and can be combined to reproduce Eq. (12.441) . In this regime, Eqs. (I2.42p and (I2.43P 
no longer represent constants of the motion. 
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3 Single-term Systems 

In special cases, a single resonant term may demonstrate a system's important dynamical at- 
tributes. This situation occurs when one resonant object is much larger than the other, and the 
larger object is on a circular orbit. In this case, the more massive resonant object remains unper- 
turbed from its orbit, and the only nonzero resonant argument contains only one nonzero value 
among the following coefficients: jfy. 

As a check on our mod el and a demonstration of its output, we reproduce some results from 



Winter and Murrayl (119971 ) ? s study of Jovian asteroid motion. Under the guises of the planar, 
circular restricted three-body problem, their study integrates Hill's equations of motion for a Jovian- 
mass planet at 1 AU and a massless asteroid at a range of values around the 2:1 commensur ability. 
The mass ratio of the Jovian planet and the central object is taken to be 10~ 3 , and the motion 
of an asteroid is computed for various and values chosen in order to demonstrate orbital 
behavior near the 2:1 commensurability. The authors fixed the origin at the system's barycenter, 
and integrated the equations of motion in a uniformly rotating frame. They also applied the common 
approximation of setting fd, and hence, to constant values. Our code does not rely on this 
approximation, which has an observable effect, even in the restricted case, on the resulting orbital 
evolution profiles. 

Figure [U illustrates an asteroid's motion with initial ci2 = 0.602 AU, ei = 0.004, A? = 0° and FIG. 1 
vp-). = 0°, with lib = {2, —1, 0, —1, 0, 0}. Direct comparison of our plots with [winter and Murray 



(Il9971 )'s Figure 14 reveal nearly equivalent profiles for the evolution of orbital elements. Differ- 
ences in W2 evolution and the presence of small ci2 modulation may be attributed to the different 
coordinate systems in which the asteroids evolved. Varying the asteroid's semimajor axis by an 
amount comparable to the barycentric correction indeed induces drastically different behavior. As 
to be expected in resonant systems, evolutionary behavior may be highly sensitive to the initial 
semimajor axes values. 

As evidenced both by observations of the Solar System and by the d'Alembert relations, one is 
less likely to come across a dynamical system in a purely inclination-based resonance where just 
one resonant angle is librating. Such a resonance must be of at least second-order, and eccentricity 
cannot play a role in the resonant evolution. Figure [2] illustrates an asteroid's motion with initial FIG. 2 
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a 2 = 0.630 AU, J 2 = 10°, A 2 = 0° and tt 2 = 60°, with lib = {4, -2, 0, 0, 0, -2} evolved over the 
same time as the asteroid in Fig. [TJ Because this resonance is second-order, the so-called "strength" 
of the resonance scales as ~ I 2 , as opposed to ~ e for a first-order eccentricity resonance. 

Figure [3] plots the relative error in the N 2 constants (Eq. 12.421) throughout the runs of the FIG. 3 
simulations in Figs. [1] and [2] for different user-inputted integration accuracy parameters of our code. 
For the eccentric asteroid in Fig. [TJ the solid line results from an accuracy parameter three orders 
of magnitude smaller than that from the dashed line. Similarly, for the inclined asteroid of Fig. 
[2J the dot-dashed and dot-dot-dot-dashed lines in Fig. [3] result from accuracy parameters which 
differ by three orders of magnitude. Figure [3] corroborates the validity of the expression for N 2 , 
and exhibits variation according only to computer precision. The relative error in the Ni constant, 
heliocentric Hamiltonian (Eq. [2722]), Jacobi Hamiltonian (Eq. I2.27p . heliocentric N 3 constant (Eq. 
12.441) . Jacobi N 3 constant (Eq. 12.451) . and angular momentum constants (Eqs. I2.19ti2.2i]) are all 
either indistinguishable from zero or are not applicable to this specific problem. 

The sole resonant argument in Figs. [1] and [2] librates about 0° and 180°, respectively. One may 
attempt to model the analytical form of the libration of resonant angles by considering Eq. (12.161) . 
Because D^ ,u ' and vary with time, and because of the summation, the equation is not easily 
integrated, and the evolution of the resona nt arg uments are sometimes by no means simple sinusoids 



(see, e.g. Fig. 17a of IWinter and Murray 



19971 ). However, I find that, in some cases, the form of 



libration profiles can be Fourier decomposed into the dominant frequencies of the system such that 

3 

4>(t) = K + J2 [Ksi-z cos (K 3i t) + K 3i _t sin (K 3i t)\ , (3.1) 

i=l 

where the summation can be taken to infinity. 

For example, the dominant frequency of the asteroid in Figs. [1] and [2] correspond to periods 
of 7.1 yr and 78 yr, as can be evidenced by considering the K 3 values from Table [A] This table TABLE 1 
provides a list of the coefficients, and their standard deviations, obtained by fitting Eq. ( 13. ip to 
the librati on profiles from Fi gs. [1] and [2j I performed the fits by using the Levenberg-Marquardt 



algorithm (IPress et al 



19921 ). We note also that the tiny standard deviations associated with the 
K 3 values indicate excellent agreement with the evolution from our code. As more Fourier terms 
are included, the fit becomes better, although the number of terms needed in some cases may be 
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prohibitive. The fit to the libration profile for the inclined asteroid is significantly better than that 
for the eccentric asteroid when comparing the ranges (2.63° vs. 26.02°) and standard deviations 
(0.14° vs. 8.16°) of the residuals. The reason for the discrepancy has to do with the different shapes 
of the libration profiles. 

The shape of some libration profiles suggest a better equation to fit with fewer terms. Such 
profiles often resemble those seen in radar Doppler velocity curves for exoplanet searches. Further, 
the computation of the evolution of the true anomaly from the eccentric anomaly (see Eq. 16. lip 
suggests that one can fit a libration profile to the following equation: 

(pit) = k + ki cos [k 2 + k 3 arctan (fc 4 tan (k 5 t))]. (3.2) 

Table |X] provides values for this equation's fit to the libration profiles in Figs. [1] and [2j The fit TABLE 2 
to the libration profile for the eccentric asteroid is a marked (one order of magnitude) improvement 
over the Fourier fit. Figure H] displays the residuals of the new fit, and the residuals have a range of FIG. 4 
3.87° and a standard deviation of 0.77°. 

4 Multiple-term Systems 

A disturbing function composed of just one resonant argument is often too simple a model for 
a realistic approximation to the evolution of bodies in a known resonance. Including the relevant 
secular terms may improve the approximation significantly, even though with added terms the 
constants of motion derived in Eqs. ( 12.421) - ( |2. 45ft strictly no longer exist. Orbit-orbit resonances 
typically involve disturbing functions with more than one resonant argument. In this case, when 
no secular arguments are considered, the Hamiltonian can be expressed as an appropriate linear 
combination of either Bi and G2, or G 3 and O4. Hence, the other four actions are constants and can 
be combined to reproduce Eqs. (12.441) or (12.451) . In this regime, Eqs. (12.421) and (I2.43P no longer 
represent constants of the motion. In the more general case, where multiple resonant and secular 
arguments dominate the motion, the Hamiltonian and total angular momentum of the system are 
the only constants guaranteed to exist. 

We expand our study of Jovian asteroid dynamics by including additional resonant and sec- 
ular arguments, thus "perturbing" our single-term models. For the eccentric asteroid of Fig. [1] 
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in the first-order 2:1 resonance, we now include, in Fig. \5[ all terms up to second-order associ- FIG. 5 
ated with the {2, —1, 0, —1, 0, 0}, {4, —2, 0, —2, 0, 0} and N arguments (dashed line), and all terms 
up to fourth-order associated with the {2,-1,0,-1,0,0}, {4,-2,0,-2,0,0}, {6,-3,0,-3,0,0}, 
{8, —4, 0, —4, 0, 0}, and K arguments (dot-dashed line). Similarly, for the inclined asteroid of Fig. [2] 
in the second-order 4:2 resonance, we now include, in Fig. [6l all terms up to second-order associated FIG. 6 
with the {4, —2, 0, 0, 0, —2} and K arguments (dashed line), and all terms up to fourth-order asso- 
ciated with the {4,-2,0,0,0,-2}, {8,-4,0,0,0,-4}, and {0,0,0,0,0,0} arguments (dot-dashed 
line). The d'Alembert relations preclude the existence of the {6, —3, 0, 0, 0, —3} argument. 

These figures demonstrate that the eccentricity and inclination profiles change significantly de- 
pending on the order of the terms included, even for the circular restricted three-body problem. 



Verasl (120071 ) demonstrates the additional consequences of eliminating the assumption that Jupiter 
remains on a circular orbit. Doing so places doubt on the validity of approximating the system by 
a single or a few terms. 

Multiple-planet systems said to be in resonance often feature one or more librating angles. The 
smaller the amplitude of these librating angles, the "deeper" into resonance the system is purported 
to be. Our model can make quantitative statements about how the depth of resonance varies with 
orbital parameters, and how these variations differ across different types of resonances. Here we 
sample phase space for regimes which could be applicable to two-planet extrasolar planetary systems. 
As a three-body problem, these systems are not likely to be restricted in many ways. However, one 
such restriction which we assume is that planets in these systems are coplanar. 

Given the extensive investigations of the 2:1 resonance, and the investigation performed in Section 

[5] of this work, here we sample phase space for a couple other relevant resonances: the first-order 

3:2, and third-order 4:1 resonances. For each, we have discovered configurations in the relatively 

small region of phase space which allows for libration of multiple arguments. The inner planet has 

a mass of 0.3Mj up (^ Mass of Saturn) and is set at 3 AU away from a Solar-mass star. The outer 

planet in the 3:2 resonant system has the same Saturn-like mass, while the outer planet in the 4:1 

system is given a mass of 0.03Mj np . Table |A] displays the "nominal" orbital configurations from TABLE 3 

which variations result in Figs. [7] -[101 The vastness of phase space allow us to sample just a small FIG. 7 

FIG. 8 

slice, presented below, and force us to defer further analysis to future studies. FIG. 9 

FIG. 10 
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These figures plot "libration amplitude," defined as half of the range of a resonant angle taken 
over a few periods. We must state such a definition because the actual evolution of these angles is 
sometimes heavily modulated and departs from a sinusoid, as in the dot-dashed profile from Fig. [51 
Angles are said to circulate if they span 360°, and do so for the vast majority of phase space. Each 
large plot symbol represents the outcome of a different simulation, and the curves which are fit 
through these symbols hence represent just an approximation to the shapes of these profiles. Figure 
[TJplots the outcomes for the 3:2 resonance, while Figs. [SlfTUlplot the outcomes for the 4:1 resonance. 
For the 3:2 resonance, asterisks and diamonds denote, respectively, the libration amplitude of the 
{3, —2, — 1, 0, 0, 0} and {3, —2, 0, —1, 0, 0} arguments. For the 4:1 resonance, crosses, squares, trian- 
gles, diamonds and asterisks, denote, respectively, the libration amplitude of the {0,0, 1, —1,0,0}, 
{4, -1, 0, -3, 0, 0}, {4, -1, -1, -2, 0, 0} {4, -1, -2, -1, 0, 0}, and {4, -1, -3, 0, 0, 0} arguments. 

Dashed lines in Figs. [7J1B] represent the "nominal" resonant semimajor axes. For both resonances, 
this value proves to be an excellent predictor for maximum depth for nearly all resonance angles. 
Note, however, the order-of-magnitude difference in a values plotted in the two figures due to the 
difference in order of the resonance. The effective "libration centers," which represent the mean 
value of libration, hover around 180° and 0°, respectively, for the asterisks and diamonds in Fig. [7J 
For Figs. [SlfTUl the libration centers are 180° for the asterisks, crosses and triangles, and 0° for the 
diamonds and squares. We don't present plots for the variation of A and w for the 3:2 resonance 
because these variations produce differences in libration widths of less than a degree. 

One can glean other results from the figures. A comparison the vertical axes reveal that in no 
case does a 3:2 resonant angle librate with an amplitude of less than 100°, in stark contrast to the 
4:1 resonance. This phenomenon most likely results from the difference in mass of the outer planet 
in the two systems; the more massive both resonant objects are, the shallower the resonance. Also, 
Figs. IHlfTOl suggest that multiple resonant arguments librate simultaneously for most of the phase 
space sampled, although occasionally some of the arguments circulate. Further, the curves drawn 
through the asterisks indicate that the libration width profiles can very drastically for different 
librating angles in the same system. Figure [10] demonstrates that in order to produce similar 
libration amplitude variations, vji had to be varied by a factor of three greater than w± did. 
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5 Dominant GJ 876 Terms 



We now consider the GJ 876 system with the outer (GJ 876 b) and inner (GJ 876 c) giant 
planets, which are thought to reside in a 2:1 resonanc e. Ever improv i ng orb i tal fits for initial system 



conditions a b ound; different sets may be found in 



Fischer et al. 



Rivera et al 



(120031 ) 



Laughlin et al. 



fl2005h . 



(120051 ) and the on-line Extrasolar Planets Encyclopedist! We chose m = O.32m , 
m 1 = 1.929m Jup , m 2 = 0M7m Jup , a x = 0.20783 AU, a 2 = 0.13031 AU, e l = 0.0251, e 2 = 0.232, 
Ai = 351.1°, A2 = 147.1°, w\ = 176.8°, and W2 = 198.3° for the numerical integration. The small 



circularized non-resonant planet G„ 

r 

evolution of the other two planets l 



876 d, at a ~ 0.02 AU, does not significantly affect the resonant 
\ We evolve the resonant planets on coplanar orbits, and thereby 
neglect the inclinations and longitude of ascending nodes. 

Figure [TT1 illustrates how the resonant planets' orbital parameters evolve using the full N-body FIG. 11 
integrator, HNbody (Rauch and Hamilton 2007, in preparation). We define 0i as the resonant angle 
corresponding to the constants {2,-1,0,-1,0,0}, and 02 as the resonant angle corresponding to 
{2,-1,-1,0,0,0}. Figure 1 demonstrates that both 0i and 02 librate about 0°, and w\ and w 2 
circulate, all with a period of ~ 9 years. The eccentricities exhibit a periodicity commensurate with 
that of the orbital angles' amplitudes of libration, whereas the semimajor axes exhibit little (1% 
- 2%) variation. The influence of the short-period terms are seen through the small modulations 
(or the noise) especially apparent in the orbital elements of GJ 876 c because of its relatively large 
eccentricity. The planets are said to be "in resonance" because of the simultaneous libration of 0i 
and 02, and GJ 876 c is said to be in "deep resonance" because of the small (< 10°) amplitude of 



A full treatment of the GJ 876 system, even to first-order in eccentricities, requires multiple 
resonant and secular arguments. Table \M lists all relevant resonant and secular terms for this TABLE 4 
system up to fourth order, and Figs. [T2iri5l illustrate the planetary evolution caused by all terms 
up to order one, two, three and four respectively. Differences from the true evolution exhibited by 
Fig. E] may be attributed to the neglect of short-period terms and even higher-order resonant and 
secular terms. The simulations in Figs. [1211151 were initialized with mean orbital element values, 



4 at |http://vo.obspnT fr/exoplanetes / encyclo / catalog. php| 

5 We have performed a full N-body simulation with HNBody (Rauch and Hamilton 2007, in preparation) which demonstrates that the 
presence of planet GJ 876 d indeed has a negligible effect. 
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obtained from averaging over the orbital evolution of a few periods of the outer planet from the 
exact numerical integration. The mean initial values used were a\ = 0.20908 AU, a 2 = 0.13008 AU, 
ei = 0.03054, e 2 = 0.2304, w x = 160.6°, w 2 = 178.6°, A x = 179.2°, and A 2 = 180.3°. 

Figure fT2l demonstrates that a full first-order treatment, which includes the leading terms from FIG. 12 
the {2, —1, —1, 0, 0, 0}, {2, —1, 0, —1, 0, 0}, {0, 0, 1, —1, 0, 0} and K arguments, forces 02 to circulate 
rather than librate and W\ to librate about 180° rather than circulate. These features are at odds 
with the full N-body integration. Therefore, analytical manipulations of these terms alone have 
limited utility when one attempts to analyze particular aspects of the system's evolution. 

The combination of terms comprising the second-order solution showcased in Fig. [13] is only FIG. 13 
marginally stable, and only for a high integration accuracy parameter (< 10~ 7 ). Nevertheless, 
amidst the chaos, both 0! and 02 do librate about 0°, but on a longer timescale (note difference 
in x axes). Short-term evolution, however, showcases fictitious asymmetric libration of varying 
amplitude about ~ 70° or ~ —70°. The switch between both libration centers occurs at extrema of 
the semimajor axes and eccentricities, and the presence of asymmetric libration centers demonstrates 



that the first-order solution is more accura te than t 
be found in the res tricted 3-body problem jBeauge 



different. Further, 



re se cond-order solution. Similar results may 



19941 ) even though that problem is inherently 



Lemaitre and Henrardl (119881 ) discuss the incorrect qualitative behavior which 



arises from improper truncation of the disturbing function as applied to asteroid studies. 

Higher-order treatments, however, not only are stable but reproduce expected features of the 

system. The third-order solution in Fig. [TH and the fourth-order solution in Fig. [15] reproduce, to FIG. 14 

FIG. 15 

varying extents, the orbital evolution shown from the N-body simulation of Fig. [TH However, in Fig. 



HI ex represents an unmodulated sinusoid, 0i lacks any modulated envelope, and 2 acquires a slight 
sawtooth form. These minor features all differ from the N-body integration, but are all partially 
remedied by the fourth-order treatment in Fig. [151 However, this same fourth-order treatment 
demonstrates a decrease in precision by altering the amplitudes and elongating the libration period 
by over a factor of two. 



The above considerations suggest that the so-called "Laplacian expansion" (described in 



Ellis and Murray 



20001 ) expansion of the disturbing function best reproduces the evolution of GJ 876 for an eccen- 



tricity order between 3 and 4. The lack of precision can be explained by the failure to satisfy the 
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Sundman criterion (ISundmanl . Il916f) for the coefficients of this expansion. In the plana r case, the 
criterion is typically expressed as (jFerraz-Melld . 1 1994 ISidlichovsky and Nesvornyl . Il994l ) : 



a 2 F(e 2 ) < ai/(ei), (5.1) 

where 



F(q) = yl + ? 2 cosh z + q + sinh z, 

(5-2) 

f(q) = V 1 + q 2 cosh z — q — sinh z, 

such that w is implicitly defined as z = gcoshz. Applying the criterion to representative GJ 876 
orbital parameters yields Fig. [T6], which demonstrates that the classic expansion will show signs of FIG. 16 
divergence at particular points in the evolution. As more resonant terms are included, the effect 
of this divergence might be amplified. Hence, the optimal order with which to model GJ 876 with 
this expansion is finite. Taken as a sequence, Figs. [T2]fT5l demonstrate that perhaps "order" is an 
inappropriate metric for classifying the accuracy of a model for systems with a phase space similar 
to that of GJ 876. We highlight these perturbative issues in order to caution future investigators 
regarding the validity of using particular resonant arguments for analytical considerations. 

Resonant systems may be affected, or even dominated, by perturbations external to the detailed 
dynamical interplay described so far. The next three sections will present potentials for different 
types of perturbations and estimate their effect on GJ 876. Table El summarizes these effects and TABLE 5 
what orbital parameters they directly affect depending on the assumptions used. 



6 Central Body Oblateness 
6.1 Overview 



The contribution of oblateness from a centra 
crucia l role in ring and small satellite dynamics. 



body to a resonant system has sho wn t o play a 



Gozdziewski and Maciejewskil ( 119981 ) and 



Shinkin 



(120011 ) h ave crea t ed an alytical models for the dynamical evolution of satellites around oblate plan- 



ets, and IWiesell ( Il982l ) has created a model of rings in resonance with an oblate planet. In this 
section, we incorporate the oblateness effect through both an averaged (over short-period terms) 
and unaveraged oblateness potential. Doing so helps determine specifically the error incurred when 
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using the averaging approximation, and helps clarify the seemingly contradictory expressions found 
in the literature. The oblateness potential, TZ^i bl \ has the form, 



n 



(obi) 



Ji I — I Pi (sin b k 



r k — V r k 



(6.1) 



where R c is the equatorial radius of the central object, Ji are the oblateness coefficients, are 
egendre polynomials, and S i s the declination with respect to the central object's equatorial plane 



Brouwer and ClemenceL 



1961 



p. 563). In order to maintain consistency with the other equations in 
this work, we wish to express Eq. (16.11) in terms of contact elements. However, we may instead use 
expressions for osculating elements because every term in the expansion for the potential depends 
only upon the resonant body positions, and not velocities (expressions for positions are identical in 
both osculating and contact elements). 



6.2 Averaged expressions 



Various different expressions for the orbital average of the oblatene ss potentia l 



rave appeared in 



Greenbergl (1198 ll ) attempted 



the literature due to differences in meaning of the orbital elements, 
to eliminate the confusion regarding these seemingly dissimilar formulas. He wrote the following 
expressions for apsidal precession rates in osculating (subscript "o"), sidereal (subscript "s"), and 
epicyclic (subscript "e") coordinates: 



w w n r 



w ~ n s 



3 fR\ 63 j2 / R\ _15 7 ( R 



2 \a s 



R 



15 



;j - - —Ja - 



R 



igM 



3 j( r Y 15 / ( R 



(6.2) 



(6.3) 



(6.4) 



Any one of 



Equation f)6.2p (IBrouwer 



:qs. (E2D-(i 



\ may be correct depending on the context of the system studied. 



19591 ) may be used wh e n the orbital elements consid e red ar e all osculating 



Keplerian elements. Equation (I6.3P (iBrouwerl . 



1946 



Murray and Dermott 



19991 . p. 270) may 



be used when dealing with observed values for the sidereal mean motion and the mean distance 
from the central object (also called the "observed semimajor axis," "apparent semimajor axis," 
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or "geometric semimajor axis"). Equation (|6.4D (lEUiot et al 



Borderies-Rappaport and Longaretti 



1994 



Murray and Dermott 



1981 



1999 



Elliot and Nicholson 



1984 



p. 269) may be used when 



the only measured orbital element is the distance from the central object, a typical characteristic 
of ring systems. 

Our model does not utilize truncated expressions for apsidal precession. Rather, we wish to 
obtain an averaged disturbing function in terms of osculating elements. Unfortunately, differing ex- 
pressions for averaged potentials have also appeared in the literature. Therefore, understanding the 
manner in which Eqs. (I6.2l) - (l6.4p are derived from disturbing functions is important in determining 
th e correct form of the a v erage d potential. 



Murray and Dermott 



(1999 



, p. 269) and 



Metria (119911 ) write expressions for averaged po- 



tentials which contain a Jf term. As the potential is a linear combination of Jj terms, aver- 
aging only over the "fast" angle known as the true anomaly would not produce a product of 
the oblateness terms from which a Jf term could be derived. Typically, Jf terms appear in 
the derivations of the variation of the orbital elements. Such derivations include binomial ex- 



pansions ( 



Kozai 



Brouwer 



1946 



Greenberg 



198 ll ). Poincare -von Z e ipel t ransformations (iBrouwer 



1959 



1962 ). Taylor expansions of orbital elements 



sions (IBorderies-Rappaport and Longaretti 



1994D 



(IKozail . 



1959), or Poincare-Lindstedt expan- 



Metrisl (ll99ll ) uses the Hori-Deprit method in 



order to illustrate how a Jf term may appear in the potential if an additional, different averaging in 
the small parameters Jj is performed over the often-neglected short-period terms. This p ro cedure 



was applied to the osculating elements in orde r to yield the mean osculatin g ones. 



resulting expression for Jf differ s from that of iMurrav and Dermott 



Further, the expression from 



Murray and Dermott 



(11999 



Metric s f|l99ll ) 



(119991 . p. 269). 



p. 269) does not contain J 3 terms 



(which were argued to be typically negligible), nor terms which are independent of eccentricity or 
inclination. These latter terms are not constant, but rather depend on the variable semimajor axis. 
The terms containing J3 vanish only if the disturbing function is averaged twice - once over the 
true anomaly, and once over the argument of pericenter. However, the latter orbital angle does not 
vary "quickly" in general. 



Given the above considerations, we include 



Metrid ' (1199ll ) Jf short-period terms, and use 



Kozai 's 



6 The elimination of the j| term in Eq. 116.41 1 is a chance cancellation; the Appendix of Borderies-Rappaport and Longaretti (1994) 
shows that the J2J4, and jS terms, for example, do not cancel. 
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(119591 ) expression for the other terms in order to obtain an averaged potential. The resulting 
expression, when truncated to "fourth-order" in eccentricities and inclinations (as is the disturbing 
function), reads: 



where, 



2 r 



i=l 



W , U) 2 , \^ ) 4 i (J ■ 2 t i W • 4 r 



+ 7j- 5) e^ sin 2 I k + T^ef, sin 2 I k cos (2ro fc - 2fl&) J 
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(6.9) 



Incorporating disturbing function partial derivatives with Eqs. (12.31) provides explicit expressions 
for the derivatives of the orbital elements for the oblateness contributions. In order to obtain the 
oblateness contribution to 0, one can combine the derivatives of the orbital elements with Eq. (12.21) . 

6.3 Complete expressions 

In the unaveraged incorporation of oblateness into our model, we use Eq. (16. ip explicitly. The 
eccentric anomaly, E k , is related to the mean anomaly through Kepler's equation: 



E k - e k sin E k = M k = X k 



w k 



(6.10) 



such that E k = E k (X k ,w k ). 
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The true anomaly is related to the eccentric anomaly through (Murray & Dermott 1999, p. 33) 



fk = fk(e k , A fe , w k ) = 2 arctan 



1 + e k 
1 - e k 



tan 



(6.11) 



Let: 



Y k = sin 4 sin (/ fc + w k - Q k ). (6.12) 

Substituting Eq. ( 16. 12[) into Eq. ( 16.11) and then incorporating the appropriate partial derivatives 
into Eq. (I2.3P allows us to determine the oblateness contribution for each orbital parameter. In 
order to compute these derivatives, the eccentric anomalies must be computed through Kepler's 
Equation numerically at each timestep. 

6.4 Application to GJ 876 

Saturn's oblateness dominates the precession of the pericenters of the system's close small satel- 



lites. 



(Ilorio 



he S un's oblateness, however, despite its significance in general relativistic computations 



20051 ). is thought to have a negligible effect on the evolution of the planets in the Solar 
System. Extrasolar systems, however, display a variety of orbital configurations which might admit 
the possibility of the central star's oblateness playing a role in the evolution. Massive (> Ms a tum) 
planets are known to orbit their stars at distances an order of magnitude shorter than the Sun- 
Mercury separation, and can evolve in orbit-orbit resonances absent from the inner Solar System. 
We further analyze the resonant planets in GJ 876 by attempting to determine under what condi- 
tions can oblateness play a role, and how likely their evolution is affected by the property of the 
star. As the effect of oblateness through inspection is typically negligible, we track the contribution 
by explicitly plotting the contribution for the rate of change of relevant orbital elements at each 
timestep. 

GJ 876 is an 0.327? ?,^ M4 dwarf with a rotational velocity sim ilar to that of the Sun (k, 2 km/s) 



(jDelfosse et al 



19981 ) and a radius just a fraction of the Sun's (jChabrier and Baraffe 



1993). The 



star's oblateness is unkno wn, likely cha n ges w ith time, and is a detailed function of the rotation 



profile, mass, and radius. 



Rozelot et al. 



measurements from several references, and 



(200 1J) compares quadrupole 



(J 2 ) and octopole (J 4 ) Solar 



Rozelot and Roeschl (119971 ) provide an upper bound to 



the Solar oblateness. As the range of values estimated for the Sun span two orders of magnitude 
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(~ 10~ 8 — 10~ 6 ), and J 2 ~ ^qRq/^q (IGodier and Rozelotl . |2000| ). where fi is the rotational 
velocity of the Sun, GJ 876's oblateness is likely to fall within this range. 

Through simulations including the 11 arguments needed for a third-order treatment, we conclude 
that the GJ 876 planets are negligibly affected by stellar oblateness. Figure [IT] plots the fractional FIG. 17 
contribution of oblateness at each timestep of w\ (solid line) and W2 (dotted line) in the unaveraged 
planar problem given J 2 = — 10~ 6 , J4 = 10 -12 and R & = 6 x 10 5 km; Fig. [TS] is equivalent FIG. 18 
except in the averaged problem. Both figures indicate an upper bound for the fractional oblateness 
contribution is 10 -6 — 10 -7 . The librational periodicity is more evident in Fig. [18] than Fig. [17] due 
to the averaged nature of the system simulated there. The magnitude of the fractional contributions 
from Figs [T7J and [18] demonstrate that any planet nonnegligibly affected by stellar oblateness must 
be much closer to the central star, at a radius which is likely to be inside the tidal circularization 
limit. 

If GJ 876 was a fast rotator (30 — 50 km/s) like other observed M dwarfs (Delfosse et al. 1998), 
then the effect of oblateness would be at least two orders of magnitude greater, and likely affect 
the long-term, if not short-term, evolution of the system. As several tens of extrasolar planets have 
been discovered within 0.1 AU of their parent stars, including several multiple systems, prospects for 
finding an oblate star with a tight resonant system are promising. Such stars would be Jupiter-like 
due to their mass, radius, and rotational period of ~ 10 hr. 



7 Central-body Precession 
7.1 Introduction 

An object's precession about its axis exerts a gravitational influence on orbiting masses. Although 



much precession-based research 



1996 



Christou and Murray 



1997 



ocuses on 



Bills 



Vlars' chaotic obliquity (e .g. 



1999 



Bouquillon and Souchay 



Hilton 



ally general models applicable to different systems are developed (e.g. 



1991 



Groten et al. 



1999 ) and although occasion 



Blitzcr 



1984 



investigators of resonant systems do not always consider the effect of precession. 



Penna 



Kozail ((I960 



999), 



conducted one of the first studie s of ho w a satellite's orbital elements are affected by the precession 



of its parent planet. 



Kinoshital (119931 ) studied th e moti on of a satellite relative to the equatorial 



plane of its oblate parent parent, and 



Rubincaml (120001 ) discusses the possibility that Pluto may 



be in "precession-orbit" resonance. Expressions for the precession contribution to the disturbing 
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functi on ar e given imp l icitly b y iGoldreichl (Il965bl ) and explicitly by iKopall (119691 ). iBrumberg et al 
(ll970l ) and lEfroimskyl (l2005bl ); using their results, we can write: 



ore) 



fik^k (1 — e|) (u>i sin I k sin Q k — uj 2 sin I k cos Q k + ^3 cos I k 



(7-1) 



where Ui, u 2 , and are the components of the precessional frequency corresponding to increasing 
moments of inertia of the central body. 

7.2 Application to GJ 876 

As displayed in Table 2, a central object's precession does not alter the semimajor axes nor 
eccentricities of the planets directly, but rather only indirectly through <fi. For GJ 876, we assume 
representative values of u)\ = uj 2 = uj 3 = 10 pHz and R Q = 6 x 10 5 km. These to values correspond 
to a precessional period of ~ 2 x 10 4 yr, which corresponds to a frequency significantly less than the 
mea n motion of the resonant bodies, a realistic regime at least in the context of precessing stellar 



jets ( INamouni 



20051 ). Figure dni plots the fractional contribution of precession at each timestep of k\ FIG. 19 



(solid line), e 2 (dotted line), w\ (dashed line), and w 2 (dash-dot line). The effect on the orbital angles 
ranges from one part in 10 -4 to one part in 10 -2 , suggesting that for planar extrasolar planetary 
systems, precession may play a greater role than central-body oblateness. Direct comparison with 
Figs. [TTirTBI demonstrates the effect of precession on the evolution of w\ and w 2 is roughly two 
orders of magnitude higher than that from oblateness. 

8 A Surrounding Thin Disk 
8.1 Introduction 

The birthplace of planets and satellites, a protoplanetary or protosatellite disk ultimately de- 
termines the dynamical interactions which occur in the system. Although capture into resonance 
is likely to and often thought to occur after disk dissipation, and requires a dissipative medium, 
captur e may be achieved while the d isk is still present, for both protoplanetary and protosatellite 



disks. 



Thommes and Lissauerl (120031 ) find that d is 



capture of planets in resonance, and iKley et al. 



A3 ca n play a crucial role in the migration and 



(120041 ) performed numerical simulations of two 



planets embedded in a thin disk, and finds resonant capture to be a common occurrence. 
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In order to approximate the mass density profile of the nascent Solar Syste m nebular disk 



a model known as t he Minimum Mass Solar Nebula (MMSN) has been devised (jWeidenschilling 



1977 



Hayashil . 



19811 ). The MMSN represents the pair of values (E , w) which, when inserted into the 
density profile, X = So(ro/r)"', provides the minimum nebular mass necessary to explain the current 
masses and positions of the planets. A variety of assumptions go into this calculation, including 
the fraction and location of nebular material that is eventually scatte r ed ou t of the system and the 



possible presence and location of the snow line (Sasseloy and Lecar 



MMSN prescription sets w = 1.5 (jWeidenschillingi . 



1977 



Hayashi 



198 lh . L 



2000). The commonly used 



We consider only the gravity of thin steady disks with a power-law density profile S = Ho(ro/r) w , 
where S and r are constant, arbitrary parameters; we take the latter to be 1 AU. Such profiles 
could represent protoplanetary or protosatellite disks. We also assume that the disk harbors two 
massive non-central bodies, each having already carved out a gap, as illustrated in Fig. [20J The FIG. 20 
relative sizes of the mass-filled regions labeled I, II, and III largely determine the dynamical evolution 
of the system. Depending on the orbital properties of planets in the system, and their interaction 
with the gas, regions I and II may be small or unoccupied with gas. We assume the disk contains 
the same total mass (= m n ) as a particular fraction, x, of the central mass, where 
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8.2 Horizontal Contributions 



w 



(8.2) 



BallabbJ (1l973l ) categorized expressions for gravitational potentials of circularly symmetric distri- 



butions of matter, in the form of ho mogeneous an d heterogeneous disks, as polynomials in semimajor 



axis. Focusing on the Solar nebula, IWardl (Il98ll ) computed averaged gravitational potentials due to 
a thin disk, and utilized Laplace coefficients in the expansions. Using the notation in the Appendix, 
we express the Laplace coefficients as (3 constants. The potential at a point distance from the 



iKu chncr (2003) extended the analysis to extrasolar systems, and derived a Minimum Mass Extrasolar Nebula (MMEN) based on 
data from 26 exoplanets in multiple-planet systems. He derived an overall steeper profile, with w = 2.0 ± 0.5, and individually fitted So 
and w values to individual multi-planet cxosystems. 
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central body in the disk plane due to the outer ("out") and inner ("in") parts of a disk, read, 
respectively, 



Pl+1 



21 V -21-w+l 



1=0 



21 + w 
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— r 



oedge J 



2/-1 r 21-W+2 



1=0 



2/-W + 2 fe 



r 
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w^2l + 2, 



.3) 



U) 



where r ie(i£ , e and r oerf9e are the edges of the entire disk, and r isap and r ogap are the boundaries of a gap 
that surrounds the point at which the potential is measured. Equations (18.31) and (18.41) represent 



just the leading-order term of the potential, not the entire potential itself. Given the schematic 
of Fig. [2H1 the disturbing functions for the outer and inner planet then become, respectively (for 
1-21 and 21 + 2), 
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for fixed values of r a , rj,, r c , r^, r e , and r/. At this stage we may choose to average over the angles 
on which T\ and r 2 are dependent; we present both approaches for the sake of completeness. Given 
the derivatives of in terms of orbital elements, disturbing function derivatives can be computed 
directly from Eqs. (I8.5l) - (l8.6p and incorporated into Lagrange's planetary equations without any 
averaging. Conversely, expressing r\ and r 2 in terms of orbital elements and averaging yields: 
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(8.7) 



where T fc and T fe are functions of and and where T k and T fc are independent of a>k,ek, 
and Jfc. The designations "oh" and "ih" represent "outer horizontal" and "inner horizontal." Gen- 
eralizing Eqs. (I8.3l) - (l8.6p to admit any value of w yields: 
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Using the binomial theorem under assumption of small e yields: 
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With the above equations, one may now determine all of the partial derivatives of the averaged 
disk disturbing functions and incorporate them into Lagrange's planetary equations. 

8.3 Vertical Contributions 
8.3.1 Overview 



Cameron and Pind (119731 ) derive ex pressions of p otential due to a cylindrical shell, but eval 



uate the expressions at the midplane. 



Wardl (1198ll ) evaluates the potential above or below the 



midplane by performing an inclination expansion and assuming that the maximum vertical (above 
the midplane) displacement of a body is much less than the distance between the disk edge and 
the body's semimajor axis (2rr's 2 <C [otfc — r '} 2 )- We adopt the more general assumption that 
2rr' s\ <C (r& — r') 2 , which is nearly equivalent except at high eccentricities. With this assumption, 
we can directly add the vertical disk contributions to the horizontal ones. In the expansion of the 
vertical disk contribution, the first term is independent of inclination and equals the corresponding 
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horizontal term exactly. The next term in the expansion is of order 0(1%) and is incorporated into 
our model. 

The vertical disk contribution is computed along similar lines as for the horizontal contribution, 
except the outer and inner potentials now read, 

Kt ut) = -W^sl (1 - cos 2 f k ) £ fj^ rf' [r^r 4-1 - r^ 1 ] , (8.14) 

1=1 

oo 0(1.5) 

= -2^E Q r »4 (1 - cos 2 f k ) £ 2l 2 w + 2 r^- 1 [r 2 ^ +2 - r 2 ^+ 2 ] , (8.15) 

i=\ 

where as previously stated, Sk = (sin Ik/2), and fk is the true anomaly. These expressions represent 
the leading order of the portion of the potential which contains inclination. Note importantly that 
the summations begin with I — 1, as opposed to those from Eqs. (I8.3l) - (l8.4p . Given the expressions 
for the disturbing functions similar to those in Eqs. (18.51) and (18.61) . and the derivatives of r&, Sk, and 
fk in terms of orbital elements, disturbing function derivatives can be incorporated into Lagrange's 
planetary equations without any averaging. Conversely, the averaged potential reads, with "ov" 
and "iv" referring to "outer vertical" and "inner vertical," 
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and, 



T fc « 4 (5 - ^4 - ^eij + ^ (44) , (8.21) 



iv (l) 1 _ 4Z+1 / 2 , l p 4\ / 2p 6 \ 

T ^ s 2 2 8 1^+2^ + Q . (8.22) 

(4 -if V(4-1) 2 V 

8.4 Application to GJ 876 

Due to the close sweeping orbits of the two resonant planets in GJ 876 and their proximity to 
the star, we consider only a disk exterior to the outer planet (r a = rj,, r c = r^, so that only the 



region labeled "III" exists). The presence of the "common gap" in between t 



re planets has been 



shown to b e the result of some numerical simulations of two planets and a disk (IBryden et al. 



2000 



KleyL 



20001 ) . An exterior disk could represent the equilibrium state remnant from the last instance 
of giant planet migration. Because an accurate treatment of the short-period variations of the 
planetary semimajor axes and eccentricities would entail the inclusion of Lindblad and corotation 
resonances for disk feedback, we focus just on the disk-induced variations of the resonant angles. 
We take r e = 0.3 AU and rj = 20 AU, but through additional simulations find that the results are 
robust against the value of r/, as little mass resides toward the outer disk edge for w > 1. The 
value of r e was chosen to lie several Hill Radii away from the apocenter of the outer planet. As 
the disk theory presented here is viable to third order, all 11 resonant and secular arguments up to 
third-order were included in the disk simulations. 

We find that disks as massive as the MMSN alters libration widths of the dominant resonant 
angles on the order of degrees and the circulation rates of the planets' longitude of pericenters on 
the order of degrees per year. The steeper the surface density profile, the greater the effect. Most 
disks perturb the planetary system enough to markedly affect the orbital angle values, but not 
enough to change their character. However, a massive or steep enough disk will warp the w profiles 
and transform libration into circulation. Figure |2~H demonstrates how this phenomenon occurs with FIG. 21 
X = 1/3, w = 3/2, r e = 0.3 AU and rf = 20 AU. The background line represents the evolution of 
02 without a disk present, and the foreground crosses were computed for the presence of a disk. Of 
further note is that the averaged disk problem fails to reproduce this circulation. The circulation is 



37 



achieved by the increase in amplitude of the short-term oscillations in the modulated envelope at 
almost 4 yrs into the orbit, and repeats with subsequent librational periods. 

As m n oc S oc TZk, and hence ultimately vj^ oc m n , a direct relationship exists between disk 
mass and the prospects for transforming libration into circulation. Equations (12.41) . (|2.12g) and 



1/2 

(I2.13ep indicate that Wk oc ttIq and Eqs. (12.91) suggest wt oc m 3 _fc. Hence, over time, as the 
disk primarily loses mass to the central star, the effect on apsidal libration is muted, and significant 
qualitative changes in apsidal behavior become more unlikely. However, transfer of disk material 
to either or both resonant bodies might maintain a significant perturbation on the apsidal angle. 
Detailed explorations of the phase space of regimes of mass loss with respect to apsidal libration 
represents a possible future avenue of study. 

9 Conclusion 

We present a model which evolves any two objects in resonance, whose orbits don't cross, using 
only the resonant and secular argument defined by the user, and, optionally, including central- 
body oblateness, central-body precession, and the presence of a surrounding nascent disk. Many 
careful studies of particular resonances rely on assumptions about what aspects of a system are 
relevant and necessary for inclusion. Our model provides a useful tool for determining quantitative 
estimates of the contributions from each argument and effect, and helps determine what must or 
should be included when analytically studying an orbit-orbit resonance. We also provide libration 
width analyses of relevant regions of phase space (Figs. l7l fT0l) . general explicit expressions for the 
variation of each resonant or secular argument considered (Eqs. I2.16H2.18]) . constants of resonant 
motion entirely in terms of orbital elements (Eqs. 12.221 12.42H2.441) . averaged oblateness potential 
terms that were previously a source of confusion in the literature (Eqs. I6.5H6.9I) . and gap-ridden 
disk potentials entirely in terms of orbital elements (Eqs. I8.7H8.13I and l8.16H8.22l) . 

We apply several aspects of this model to the GJ 876 extrasolar planetary system, and conclude 
that at least a third-order treatment is necessary in order to mimic the qualitative evolution of the 
system. A fourth-order treatment improves the approximation in some but not all respects due 
to Sundman's convergence criteria for Laplacian expansions, bringing into question the viability of 
using "order" as the metric for quantifying the accuracy of a resonant system. The GJ 876 planets 
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are negligibly influenced by the oblateness and precession of the central star, but suggest that given 
the right conditions in other exosystems, these effects may play a role in the dynamical evolution. 
A protoplanetary disk varies, sometimes significantly, librational amplitudes and circulation rates 
of extrasolar systems, and may even convert one type of motion into another. Such dynamical flags 
may constrain unknown orbital parameters in newly discovered exosystems. 
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A Appendix 



The quantity ff p) (Eq. I2TTUD is a function of both semimajor axes through the Laplace coeffi- 



cients, b { P. W 



len expressed as an infinite h yperKeo metric series, the coefficient s may be expressed 



as, for j 7^ 0, (IBrouwer and Clemence 



1961 



p. 495; 



Murray and Dermott 



1999 



p. 237): 



{\j\,s) 



(b'[.«)„|i|+2 j_ R {\3\,s) a \j\+i 



E 

v=l 



(A.l) 



such that 



0» 



2 2 ~ l (2s + 2i-4)H 



(z-l)!(|j|+z-l)!( 2s -2)H 
For the case of j = we can better compute the coefficients by using the following formula: 



v=s+i+\j\-2 

n » 



(A.2) 



n 

«=2 



s + w-2 
1 



, i > 2 



(A.3) 



with P 



(o,») 



2. 



By expressing as a polynomial in a = (12/cii, we can then take the necessary analytical partial 
derivatives of C^'^ and avoid integration when computing Laplace coefficients. The differential 
operator T> acts on a such that for yth-order eccentricity-type resonances, 



Y 



y=0 y=0 



v=l 



■y+l 

H(\j\+2v-i) 

i=2 



^'l4) Q! ]i|+2,-2 ) 



(A.4) 



where represents the function of j 
may be read off from Appendix B of 



corresponding to the p articu lar term in the resonance which 



Murray and Dermott 



(1999). The term in square brackets 



equals unity when y = 0. For each value of v , the corresponding coefficient of a may be computed. 
For inclination resonances, we also need to compute: 



T Y-l 



T Y-l 



£ £ ^f'^-V = E E ^ E 

t=l y=l i=l y=0 v=l 



v — 1) 



i=2 



where i is a counter for the linear combination of up to T terms which may appear in the expression 
for a given inclination term (for example, fg). Similarly, we have: 
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T Y-2 



T Y-2 



EE^'^-^EE^'E 



t=\ y=2 



t=l y=2 



v=l 



J/-1 



Y[(\j\ + 2v-i) 



,i=2 



We can finally express «f' p \ for a particular term, as: 



(A.6) 



ft 



o, 



2 



2 



for Z < \j\ 
for Z = \j\ 

for Z — I j I > and odd 



PHI IIS (* - i) , for Z — Ul > and even. 

2 



(A.7) 
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Figure and Table Captions 

Figure [U Evolution of a Jovian asteroid on an eccentric orbit. The resonant angle 0ub = 
{2,-1,0,-1,0,0}. 

Figure [2j Evolution of a Jovian asteroid on an inclined orbit. The resonant angle 0ub = 
{4,-2,0,0,0,-2}. 

Figure [3j Relative error in N2 for the asteroid evolution in Fig. [1] for user-inputted integration 
accuracy parameters of 10~ 8 (solid) and 10~ 5 (dashed), and for the evolution in Fig. [2] for accuracy 
parameters of 10~ 8 (dashed-dot) and 10" 5 (dashed-dot-dot-dot). 

Figure HI Deviation of the libration profile in Fig. [I] from the model of Eq. (13. 2p using coefficients 
in Table \M 

Figure [5j The eccentricity profile in Fig. [1] (solid line) and those derived with all terms up 
to second-order from the arguments {2,-1,0,-1,0,0}, {4,-2,0,-2,0,0}, and K (dashed line) 
and with all terms up to fourth-order from the arguments {2, —1, 0, —1, 0, 0}, {4, —2, 0, —2, 0, 0}, 
{6, -3, 0, -3, 0, 0}, {8, -4, 0, -4, 0, 0}, and K (dot-dashed line). 

Figure [6j The inclination profile in Fig. [2] (solid line) and those derived with all terms up to 
second-order from the arguments {4, —2, 0, 0, 0, —2} and K (dashed line) and with all terms up to 
fourth-order from the arguments {4, —2, 0, —2, 0, 0}, {8, —4, 0, —4, 0, 0}, and K (dot-dashed line). 

Figure [Tj Libration amplitudes for two planets in a 3:2 resonance with parameters that vary 
from the state given in Table [A] Diamonds and asterisks represent half the variation exhibited 
by the {3, —2, 0, —1, 0, 0} and {3, —2, —1, 0, 0, 0} arguments respectively for individual simulations. 
The dashed line indicates a "nominal" a value of (3/2)^ 2//3 \ 

Figure El Libration amplitudes for two planets in a 4:1 resonance with parameters that vary from 
the state given in Table [A] Crosses, squares, triangles, diamonds and asterisks represent half the 
variation exhibited by the {0, 0, 1, -1, 0, 0}, {4, -1, 0, -3, 0, 0}, {4, -1, -1, -2, 0, 0}, {4, -1, -2, -1, 0, 
{4,-1,-3,0,0,0} arguments respectively for individual simulations. The dashed line indicates a 
"nominal" a value of (4/l)(" 2 / 3 ). 
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Figure [HI Same as Fig. [S]but where Ai and A 2 are varied. 
Figure O Same as Fig. M but where w\ and w 2 are varied. 

Figure [TTJ GJ 876 b's and c's evolution from a full N-body integration. The angles <pi and (f) 2 
undergo libration about 0°, and w\ and w 2 undergo circulation. 

Figure [12) GJ 876 b's and c's evolution due to all resonant and secular terms up to first order 
in eccentricities. 

Figure [T3l GJ 876 b's and c's evolution due to all resonant and secular terms up to second order 
in eccentricities. 

Figure [TU GJ 876 b's and c's evolution due to all resonant and secular terms up to third order 
in eccentricities. 

Figure [151 GJ 876 b's and c's evolution due to all resonant and secular terms up to fourth order 
in eccentricities. 

Figure [TBI Laplacian convergence region for minimum values of a% and a 2 (dashed lines), maxi- 
mum values of a± and a 2 (dash-dot lines), minimum value of a\ and maximum value of a 2 (dotted 
lines), and maximum value of ai and minimum value of a 2 (solid lines), for eccentricity values (dots) 
achieved during the simulation. 

Figure [171 Central-body oblateness contribution to system evolution per timestep expressed as 
a fraction of the unperturbed values of w\ (solid line) and w 2 (dotted line) assuming values of 
J 2 = — 10~ 6 , J4 = 10~ 12 and R Q = 6 x 10 5 km in the unaveraged problem. 

Figure Central-body oblateness contribution to system evolution per timestep expressed as 
a fraction of the unperturbed values of w\ (solid line) and w 2 (dotted line) assuming values of 
J 2 = — 10~ 6 , J4 = 10~ 12 and R Q = 6 x 10 5 km in the averaged problem. 

Figure [19) Precessional contribution to system evolution per timestep expressed as a fraction of 
the unperturbed values of k\ (solid line), e 2 (dotted line), w\ (dashed line), and rij 2 (dash-dot line). 
Values of W\ = uj 2 = w 3 = 10 pHz and R Q = 6 x 10 5 km are assumed. 
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Figure [2UJ Cartoon of two planets embedded in a surrounding thin disk with gaps. 

Figure [2TJ Evolution of <p2 without a disk (background curve) and with a disk (foreground 
crosses) that is 1/3 the mass of the central star with a power law exponent of w — 1.5 and with 
r e = 0.3 AU, rf = 20 AU. Note how the disk converts resonant libration into circulation. 

Table 1: Constants and their standard deviation obtained in Eq. (13. IP 's analytical fit to the 
libration profile for the eccentric asteroid from Fig. [T] and the inclined asteroid from Fig. [2J The 
range and standard deviation of the difference between the fit and the model for the eccentric 
asteroid is 26.02° and 8.16°, and for the inclined asteroid is 2.63° and 0.14°. 

Table 2: Constants and their standard deviations obtained in Eq. (I3.2p 's analytical fit to the 
libration profile for the eccentric asteroid from Fig. [T] and the inclined asteroid from Fig. [2J The 
range and standard deviation of the difference between the fit and the model for the eccentric 
asteroid is 3.87° and 0.77°, and for the inclined asteroid is 2.95° and 0.18°. 

Table 3: Initial parameters which are varied to produce Figs. IT1U01 

Table 4: Summary of resonant and secular disturbing function cosine arguments, each of which 
may admit more than one term, relevant to the GJ 876 system up to fourth-order in eccentricities. 

Table 5: Summary of the orbital elements directly affected (*) by central-body oblateness and 
precession, and a nascent disk, under the assumptions of orbit averaging and planarity. 
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Table 1 



Constants 


Eccentric asteroid 


Eccentric asteroid a 


Inclined asteroid 


Inclined asteroid o 


K 


-0.204 


0.10 


180 


0.0053 


K x 


1.15 


0.24 


60.6 


0.0087 


K 2 


38.0 


0.14 


11.5 


0.017 


K 3 


0.883 


0.000061 


0.0805 


0.0000028 


K A 


0.870 


0.127 


-0.290 


0.015 


K 5 


0.041 


0.253 


-0.556 


0.010 


K 6 


1.60 


0.0030 


0.238 


0.00027 


K 7 


0.0494 


0.127 


-0.409 


0.0077 


K s 


-0.00379 


0.25 


0.0924 


0.016 


K 9 


3.23 


0.051 


0.406 


0.00041 
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Table 2 



Constants 


Eccentric asteroid 


Eccentric asteroid a 


Inclined asteroid 


Inclined asteroid o 




-0.0694 


0.0098 


180 


0.014 


h 


36.8 


0.013 


-61.3 


0.015 


fc 2 


-1.56 


0.00053 


2.95 


0.00039 


fc 3 


1.97 


0.00047 


2.03 


0.00086 




2.05 


0.0020 


0.964 


0.00061 




0.442 


0.0000020 


0.0402 


0.0000018 
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Table 3 



Resonance 


mi 


m 2 


a 


ei 


S2 


Ai 


A 2 




tu 2 


3:2 


0.3m Jup 


0.3mj 


0.763 


0.05 


0.15 


281 


210 


140 


340 


4:1 


OMmjup 


0.3m.j 


0.397 


0.25 


0.01 


327 


220 


186.5 


340 
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Table 4 



"resonant" (R) or "secular" (S) 


argument 


lst-order 


2nd-order 


3rd-ordcr 


4th-order 


S 







p 2 p 2 

e li e 2 




e}, e|, e 2 e 2 


S 






eie 2 




efe 2 , eie^ 


S 


2w\ — 2u72 








e l e 2 


R 


2Ai — A2 — zui = & 


ei 




e l: e l e 2 




R 


2Al - A 2 - -K72 = 4>1 


e2 




el, eje 2 




R 


4Ai - 2A 2 - 2zu 1 




4 




4 22 

e l; e l e 2 


R 


4Al — 2A2 — TUl — TU2 




eie 2 




eie|, efe 2 


R 


4Ai - 2A 2 - 2tu 2 




e 2 




4 22 

e 2i e l e 2 


R 


6A1 — 3A2 — 3tui 






e\ 




R 


6A1 — 3A2 — 2zui — W2 






efe 2 




R 


6Al — 3A2 — TDl — 2lJJ2 






e\e\ 




R 


6A1 — 3A2 — 3ZU2 






e 2 




R 


8A1 - 4A 2 - 4tui 








e l 


R 


8A1 — 4A2 — 3n7i — ZU2 








e?e 2 


R 


8A1 — 4A2 — 2tjj\ — 2w2 








e l e 2 


R 


8Al — 4A2 — H7l — 3n72 








e\e\ 


R 


8A1 - 4A 2 - 4tu 2 








e 2 
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Table 5 



effect: 




oblateness 


disk 


precession 


averaged? 




yes 


yes 


no 


no 


yes 


yes 


no 


no 


no 


planar? 




yes 


no 


yes 


no 


yes 


no 


yes 


no 


no 




d 






















e 






★ 
















j 






















e 






















w 






















n 
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Figure 2: 
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